require(ggplot2)
require(ggthemes)


rm(list=ls())


### This script will read in all the midline traces from stickleback C-start videos, find the points 1/3rd, 1/2, and 2/3rds down the body,
  # then find the Menger curvature of these three points. This is the curvature of the unique circle that fits these three points.
  # I use these points because they approximately correspond to the position of the three points in the 3-point bending experiment.
  # I measure curvature in that experiment in the same way, so these two results should be comparable.

  

dir <- "Cstart Midline curves"
setwd(dir)


#Read the names of the .csv files in the current folder
fileNames <- list.files(pattern="*.csv", all.files=TRUE, no..=TRUE)
#Open all the files
list2env(
  lapply(setNames(fileNames, make.names(gsub("*.csv$", "", fileNames))), 
         read.csv), envir = .GlobalEnv)

#Preallocate a data frame to hold all the curvature data
finalCurveData <- data.frame(Date=character(),Population=character(), Treatment=character(),Animal=integer(),Sequence=integer(),maxCurve=double())
#This loop open each file, performs the analysis, then saves out the data in a new folder

for(i in 1:length(fileNames)){
  nameNow <- fileNames[i]
  #Get rid of ".csv" in the filename
  nameNow2 <- substr(nameNow,1,nchar(nameNow)-4)
  #Add an x before the filename (because R data frame names can't start with a number)
  nameNow3 <- paste('X',nameNow2,sep='')
  #Open the file
  Data <- eval(parse(text= nameNow3))
  
  #Extract video data from filename
  splitStrings <- strsplit(nameNow2,fixed=FALSE, split = '_')
  dateNow <- splitStrings[[1]][1]
  popNow <- splitStrings[[1]][2]
  treatNow <- splitStrings[[1]][3]
  aNow <- as.integer(substr(splitStrings[[1]][4],2,3))
    TANow <- paste(treatNow,aNow,sep='')
  sNow <- as.integer(substr(splitStrings[[1]][5],2,3))
    
  #This loops takes the current video and finds the curvature of the middle third of the body in each frame.
  mCurves <- vector(length=ncol(Data)/2)
  for(k in seq(1,ncol(Data),2)){
    
    #Convert pixel coordinates to centimeters. Calibration value is 75.0 pixels/cm
    dataCal <- Data/75
    
    #Only keep coordinates for the current frame
    xNow <- dataCal[,k]
    yNow <- dataCal[,k+1]
    
    #Standardize body length to 1 by finding the total length of the body curve, centering the curve about the origin, and dividing coordinates by body length
    distances <- vector(length=length(xNow)-1)
      for (j in 1:length(xNow)-1){
        distances[j] <- sqrt((xNow[j+1]-xNow[j])^2+(yNow[j+1]-yNow[j])^2)
      }
    
    cumSumLength <- cumsum(distances)
    stLenNow <- cumSumLength[length(cumSumLength)]
    
    xNowSizeCal <- xNow/stLenNow
    yNowSizeCal <- yNow/stLenNow
    
    
    #Only keep the middle 1/3rd of the body because this is where the coelom (and fibrosis) is. Also to eliminate artifacts near the head/tail
      #just use it for plotting
    xNowCal3rd <- xNowSizeCal[c(68:133)]
    yNowCal3rd <- yNowSizeCal[c(68:133)]

    
    #plot(xNow,yNow, xlim=c(0,5),ylim=c(-2.5,2.5))
    #plot(xNowSizeCal,yNowSizeCal, xlim=c(0,5),ylim=c(-2.5,2.5))
    #plot(xNowCal3rd,yNowCal3rd, xlim=c(0,5),ylim=c(-2.5,2.5))
     

     
     #Find the Menger curvature of the triangle made by the three points 1/3rd, 1/2, and 2/3rds of the way down the body.
        #Menger curvature is the curvature of the unique circle made by any set of 3 points. Curvature is 1/radius.

       
       x1 <- xNowSizeCal[67]
       x2 <- xNowSizeCal[100]
       x3 <- xNowSizeCal[133]
       
       y1 <- yNowSizeCal[67]
       y2 <- yNowSizeCal[100]
       y3 <- yNowSizeCal[133]
       
       length1 <- sqrt((x2-x1)^2+(y2-y1)^2)
       length2 <- sqrt((x3-x2)^2+(y3-y2)^2)
       length3 <- sqrt((x1-x3)^2+(y1-y3)^2)
       
       areaTri <- abs((x1*(y2-y3)+x2*(y3-y1)+x3*(y1-y2))/2)
       
       #if area=0 then the points are co-linear and the curvature will be infinite, so we set it to NaN instead. Otherwise, we calculate curvature
       if (areaTri==0) {
         mCurveNow <- 'NaN'
       } else {
         mCurveNow <- (4*areaTri)/(length1*length2*length3)
       }
       
     mCurves[(k+1)/2] <- mCurveNow
     maxCurve <- max(mCurves)
     
     tempDF <- data.frame(Date=dateNow, Population=popNow, Treatment=treatNow,Animal=TANow,Sequence=sNow,maxCurve=maxCurve)
     
  }
  
  finalCurveData[i,] <- tempDF
  
  
}


###Find max values of the curvatures

#show max values for each individual so I can enter it into spreadsheet

allIDs <- unique(finalCurveData$Animal)
for (k in 1:length(allIDs)) {
  
  temp <- max(finalCurveData$maxCurve[finalCurveData$Animal==allIDs[k]])
  temp2 <- paste(allIDs[k],' ',temp)
  print(temp2)
}

#Put max values in data frame
#make new columns to hold only max values
naVect <- rep(NA,nrow(finalCurveData))
naDF <- data.frame(maxBodCurve=naVect)
DataWithMax <- cbind(finalCurveData,naDF)

#Find the maxima, find where they are in the data frame, and copy the value to a new column for that metric
#doesn't work for duration because it is copying all instances that match each maximum duration value.
#I don't care because duration isn't a performance metric, more of an effort metric.
for (k in 1:length(allIDs)) {
  for (j in c(6,7)){
    tempPerf <- max(DataWithMax[DataWithMax$Animal==allIDs[k],j])
    DataWithMax[DataWithMax[,j]==tempPerf,j+3] <- tempPerf
    
  }
}




#Save the dataframe
write.csv(DataWithMax, file="finalCurveData_April20.2022.csv")
